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ABSTRACT 

We investigate the formation and evolution of giant molecular clouds (GMCs) in a Milky- Way-like 
disk galaxy with a flat rotation curve. We perform a series of 3D adaptive mesh refinement (AMR) 
numerical simulations that follow both the global evolution on scales of ^ 20 kpc and resolve down to 
scales < 10 pc with a multiphase atomic interstellar medium (ISM). In this first study, we omit star 
formation and feedback, and focus on the processes of gravitational instability and cloud collisions and 
interactions. We define clouds as regions with nu > 100 cm"'^ and track the evolution of individual 
clouds as they orbit through the galaxy from their birth to their eventual destruction via merger or 
via destructive collision with another cloud. After ^ 140 Myr a large fraction of the gas in the disk 
has fragmented into clouds with masses ~ 10^ Mq and a mass spectrum similar to that of Galactic 
GMCs. The disk settles into a quasi steady state in which gravitational scattering of clouds keeps the 
disk near the threshold of global gravitational instability. The cloud collision time is found to be a 
small fraction, ~ 1/5, of the orbital time, and this is an efficient mechanism to inject turbulence into 
the clouds. This helps to keep clouds only moderately gravitationally bound, with virial parameters of 
order unity. Many other observed CMC properties, such as mass surface density, angular momentum, 
velocity dispersion, and vertical distribution, can be accounted for in this simple model with no stellar 
feedback. 

Subject headings: galaxies: spiral, galaxies: ISM, galaxies: star clusters, methods: numerical, ISM: 
structure, ISM: clouds, stars: formation 
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1. INTRODUCTION 

The formation of stars in disk galaxies is one of 
the most basic processes controlling galactic evolution. 
While there are many other important effects, such as 
galactic interactions and infall of diffuse gas, ultimately 
a large fraction of gas settles in rotationally supported 
disks, where the majority of the stellar population is 
born. The appearance of many galaxies is dominated 
by the light from massive, short-lived stars whose spatial 
distribution is controlled by this star formation process. 
The injection of heavy elements from winds and super- 
novae will occur mostly from this disk environment. 

Empirical correlations have been found between the 
star formation rate (SFR), gas content and global galac- 
tic dynamics. Based on a sample of about 100 nearby 
galaxies and circumnuclear starburst disks, Kennicutt 
(1998) found relatively simple relations between _the 
globally-averaged star formation rate per unit area, Sgfr, 
and the total (H I and H2) gas mass surface density, Egas- 
They can be related via 



Egfr — ^sfrEg"''', 



(1) 

with Asft. = (2.5 ± 0.7) X 10"'* Mq yr'^ kpc'^ and a^h = 
1.4 ± 0.15. Alternatively, an equally good fit to the data 
is given by _ 

Esfr — -BgfrEgriout J (2) 

where -Bsfr — 0.017 and rJout is the orbital angular fre- 
quency at the outer radius that is used to perform the 
disk averages. This last relation implies that a fixed frac- 
tion, about 10%, of the gas is turned into stars every 
outer orbital timescale of the star-forming disk. 

Martin & Kennicutt (2001) argued that the outer edge 
of the star-forming disk is set by the boundary between 



the gravitationally unstable Q ^ 1 inner disk and the 
gravitationally stable Q ^ 1 outer region, where Q is the 
Toomre (1964) stability parameter 



Q 



(3) 



where <jg is the ID gas velocity dispersion in the disk 
plane and k is the epicyclic frequency: 

'^'^V2^(H_/3)i/2. (4) 



1 



r dvc 
dr 



Here, Vc is the circular velocity at a particular galac- 
tocentric radius r, and (3 = d Invc/d Inr, which is 
for a flat rotation curve. The precise value of Q be- 
low which gas becomes unstable is also affected by the 
destabilizing influence of the potential due to a stellar 
disk (Jog & Solomon 1984, Jog 1996, Kim & Ostriker 
2007). Schaye (2004) argued that the edges of star form- 
ing disks may rather be set by the location where much 
of a galaxy's atomic interstellar medium transitions from 
the cold (- 300 K) to the warm (T ~ 8000 K) phase. 

When galaxies are examined on smaller (^kpc) scales, 
e.g. by taking azimuthal averages, star formation rela- 
tions similar to eqs. ((1]) & (21) are found (Wong & Blitz 
2002, Kennicutt et al. 2007, Leroy et al. 2008, Bigiel et 
al. 2008), although there is still debate as to whether it is 
the total or just the molecular gas mass that is most fun- 
damental for controlling the star formation rate (Blitz & 
Rosolowsky 2006, Leroy et al. 2008, Bigiel et al. 2008). 

1.1. GMCs: Formation, Evolution, Observed Properties 

As we move to smaller scales (~ 100 pc) we see, 
based mostly on surveys of CO emission in the Milky 
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Way (Solomon et al. 1987; Dame, Hartmann & Thad- 
deus 2001; Jackson ct al. 2006) that most of the dense, 
cold, potentially star-forming gas is organized into giant 
molecular clouds (GMCs) with masses ^ 10^ Mq and 
radii ~ 30 pc (with these particular values of mass and 
size corresponding to volume-averaged number densities 
of H nuclei of nu = 260 cm^'' for a spherical cloud). 
Given that most of the Galactic molecular gas is orga- 
nized in GMCs and that GMCs have approximately equal 
mass atomic envelopes (Blitz 1990), the results of Wolfire 
et al. (2003) for the radial distribution of Galactic molec- 
ular and atomic gas indicate that a large fraction, ^ 1/3, 
of the total gas mass in the Milky Way inside the solar 
circle is associated with these structures. Many questions 
concerning GMCs are still debated, including whether 
they are long or short-lived compared to their free fall 
timescale, 



Bertoldi & McKee (1992) define the virial parameter 
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4.35xl0«(— ^) '^'yr, (5) 
V 100 cm / 



and whether they are gravitationally bound (McKee & 
Ostriker 2007, and references therein). 

Observations of GMCs are most complete in the Milky 
Way interior to the solar circle (Solomon et al. 1987; 
Wilhams & McKee 1997; Heycr et al. 2008). Except 
for the vicinity of the Galactic center, there is a dearth 
of clouds in the inner ~ 2 kpc. GMCs occupy a very 
thin vertical distribution in the Galaxy: Bronfman et al. 
(2000) derived a half- intensity height of Zi/2 — 60 pc; 
Stark & Lee (2005) derived a vertical scale height of < 
35 pc. 

Stark & Brand (1989) derived a ID cloud to cloud RMS 

velocity dispersion of 7.8 ± 0.5 kms^-'^ from a study of 
GMCs within 3 kpc of the Sun, although this estimate 
includes small scale streaming motions. 

Based on the ^^CO survey data of Solomon et al. 
(1987), Williams & McKee (1997) derived a cloud mass 
function of the form 



d\nMr 



(6) 



for Mc < Mu, with dAfc{Mc) being the number of clouds 
with masses in the range [Mc,Mc{l + dlnMc)]. Based 
on clouds in the mass range of a few x 10^ Mq to a few 
xlO^ Mq, Wilhams & McKee (1997) estimated that the 
population of GMCs inside the solar circle is described by 
ttc = 0.6, Men = 63 and M„ = 6 x 10^ Mq. If the obser- 
vational surveys have higher degrees of incompleteness at 
lower masses, Williams & McKee estimated ac = 0.85, 
A/'cu = 25 and M„ = 6 x 10^ Mq. Note these mass es- 
timates did not include any atomic gas that might be 
associated with the GMCs, which could be a substantial 
fraction (^ 1/2) of the total mass bound to the cloud 
(Blitz 1990). 

The mean mass surface density, Sc, of Galactic GMCs 
was derived by Solomon et al. (1987) to have a me- 
dian value of about 200 Mq pc^^. Recently, Heycr ct al. 
(2008) derived a median value of 42 Mq pc^^, based on 
an LTE analysis of ^^CO survey data. They estimate 
the true values are larger because of subthermal exci- 
tation and abundance variations and arc in the range 
80 — 120 Mq pc^^. Again these estimates do not include 
any atomic gas associated with the clouds. 



GMc 



(7) 



where CTc is the mass-averaged ID velocity dispersion of 
the cloud, which includes thermal and nonthermal con- 
tributions via (Tc = (Cg + af^^ ^Y^^ where cTnt.c is the ID 
RMS velocity dispersion about the cloud's center of mass 
velocity, and Rc,a is a measure of cloud radius based on 
its projected area. The virial parameter is a measure 
of the ratio of the kinetic to gravitational energy of a 
cloud. For a imiform, spherical cloud, avir = 1 implies 
the total kinetic energy of the cloud is half the magnitude 
of the gravitational energy. Nonspherical and nonuni- 
form density distributions typically make only modest ef- 
fects, parameterized via the dimensionless factor a in the 
gravitational energy equation W — —{3/5)aGM^/Rc,A- 
Bertoldi & McKee (1992) estimate that cloud major to 
minor axis ratios of a factor of 3 only cause a to deviate 
from unity by < 10%, while a power law density profile 
p oc r~^/^ results in o = 1.25. Observationally, Heyer et 
al. (2008) find a median value of ctvir of about unity. 

Phillips (1999) considered the rotational properties of 
Galactic GMCs, finding a significant spread in the direc- 
tions of the angular momentum vectors (as derived from 
radial velocity gradients and as projected on the plane of 
the sky), indicating that a substantial fraction of GMCs 
rotate in a retrograde direction with respect to Galactic 
rotation. 

Observations of GMCs in other galaxies tend to find 
they have similar properties as Galactic GMCs, including 
their velocity dispersions and virial parameters (Bolatto 
et al. 2008). Rosolowsky et al. (2003) measured the 
rotational properties of large GMCs in M33, finding rel- 
atively small values of specific angular momentum, and 
a substantial population (~ 40%) of clouds with rotation 
directions retrograde with respect to that of the galaxy. 
Fukui et al. (2008) studied 164 GMCs in the Large Mag- 
ellanic Cloud, concluding these clouds were also close to 
virial equilibrium. 

The star formation activity within GMCs appears 
highly clustered. Most stars are born from star-forming 
clumps that turn into star clusters (Lada & Lada 2003), 
with initial radii ^ 1 pc. Here the local overall star 
formation efficiency is relatively high, e* = M*/Mgas 
0.1 — 0.5, where M^, is the total mass of stars formed from 
the mass of gas Mgas that occupies the same volume as 
forming star cluster. Conversely, most of the volume of 
GMCs is not actively forming stars, perhaps because of 
magnetic field support (Crutcher 2005). 

Theoretical work on the formation of GMCs has led 
to different schools of thought (McKee & Ostriker 2007, 
and references therein). 'Top-down' formation mecha- 
nisms suggest GMCs form via large-scale gravitational 
or magnetic disk instabilities (e.g. Kim, Ostriker, & 
Stone 2003; Shetty & Ostriker 2006; Glover & Mac Low 
2007), whereas 'bottom-up' processes have GMCs form- 
ing via agglomeration from inelastic collisions of clouds 
(e.g. Kwan 1979) or from turbulent or colliding flows 
(e.g. Bergin et al. 2004; Heitsch et al. 2008). It is possi- 
ble that both processes may be important depending on 
the galactic environment (Dobbs 2008). If GMCs are rel- 
atively long-lived and contain a large fraction of the total 
ISM, then the formation and destruction of GMCs to and 
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from the atomic phase may be less important than the in- 
teractions between those aheady existing (Tan 2000): in 
other words, the nature of a particular cloud may change 
more drastically via merging collisions with other GMCs 
than via flows into or out of the atomic phase. 

In this paper, as a first step towards understanding 
galactic star formation rates, we will examine the forma- 
tion and evolution of GMCs in the context of the global 
dynamics of a flat rotational curve disk galaxy. 

1.2. Theories of Disk Galaxy Star Formation Rates 

Numerous theories have been proposed to explain the 
observed galaxy-scale, Kennicutt star formation rela- 
tions. One group of theories is based on the growth rate 
of gravitational perturbations in a disk. The timescale 
for perturbation growth can be expressed as igrow c>c 
(G/9g)-°-5 (e.g. Larson 1988; Elmegreen 1994; Wang & 
Silk 1994), and so psfr oc Pg/igrow c>c Pg'^. Assuming a 
constant disk scale height, we obtain equation ([T]) with 
asfr = 1.5 for local mass surface densities. However, disk- 
averaged quantities will depend on the radial gas distri- 
bution. We can also express igrow ~ '^g / {''^GYig) ~ Q/k. 
Perturbation growth via swing amplification in a differ- 
entially rotating disk occurs at a similar rate (e.g. Lar- 
son 1988). By assuming star formation self-regulates and 
keeps Q constant, Larson (1988) and Wang & Silk (1994) 
predicted Sgfr oc Sg/tgrow « Sgil, since k oc fi, for 
disks with flat rotation curves. Li, Mac Low & Klessen 
(2006) presented isothermal smooth particle hydrody- 
namic (SPH) simulations of disk galaxies, from which 
they concluded rate of the nonlinear development of grav- 
itational instability determines the local and global Ken- 
nicutt relations. 

However, these theories that involve the star formation 
rate being set by the growth rate of large scale gravita- 
tional instabilities leading to CMC formation, have dif- 
ficulty in explaining why a large mass fraction of the gas 
is already organized into GMCs. If the GMCs are form- 
ing rapidly in '--^ 1 dynamical timescale the mass flux 
into (and out of) GMCs would be ~ 100 times greater 
than that from GMCs into star clusters (Zuckerman & 
Evans 1974). To understand global star formation rates, 
it seems more reasonable to look for processes that cre- 
ate the star-forming parsec-scale clumps within GMCs, 
rather than the processes that create the GMCs them- 
selves. 

The spatial correlation of star formation with large 
scale spiral structure in some disk galaxies motivates 
theories for the triggering of star formation during the 
passage of gas through density waves. Wyse (1986) and 
Wyse & Silk (1989) propose a SFR law of the form 

SsfrCxS^=''(r!-r!p) (8) 

where VLp is the pattern frequency of the spiral density 
wave. In the limit of small fip and for ash = 1 we re- 
cover equation ([2]). Increased cloud collision rates and 
increased perturbation growth rates in the arms, where 
Q is locally lowered, have been suggested as the star for- 
mation triggering mechanism (e.g. Dobbs 2008). Kim, 
Kim, & Ostriker (2008) examined the development of 
thermal instabilities in galactic spiral shocks. 

One prediction of these theories is a correlation of 
SFR with the density wave amplitude. However this is 



not observed (Elmegreen & Elmegreen 1986; McCall & 
Schmidt 1986, Kennicutt 1989). Furthermore such theo- 
ries have difficulty explaining star formation in galaxies 
where there is a lack of organized star formation features, 
as in flocculent spirals (e.g. Thornley & Mundy 1997a, 
1997b; Grosbol & Patsis 1998). GMCs are present and 
SFRs are similar to those systems where star formation is 
organized into spiral patterns. This suggests stellar disk 
instabilities, which create spiral density waves, and gas 
instabilities, which lead to GMCs and large-scale star for- 
mation, arc decoupled (Kennicutt 1989; Seiden & Schul- 
man 1990). 

Tan (2000) proposed a model in which the majority 
of star formation in disk galaxies occurs in the pressur- 
ized regions triggered by GMC-GMC collisions. GMCs 
are observed to occupy a very thin vertical distribution 
in the Galaxy (Stark & Lee 2005), which is similar to 
the actual sizes of the clouds. In this effectively two di- 
mensional system, collisions are set by galactic shear at 
impact parameters of about the tidal radius of the clouds 
(Gammie, Ostriker, & Jog 1991), and for a Q ~ 1 disk 
with a relatively large mass fraction, ~ 1/2, in GMCs 
and associated gas, collisions are expected to occur ap- 
proximately every 20% of an orbital time. In this way, 
for flat rotation curve galaxies, equation ([2]) is recovered, 
and the cloud collision mechanism is the link between 
the global galactic dynamics and the parsec-scale star- 
forming clumps of GMCs. One prediction of this model 
is a reduction in the star formation rate per orbital time 
for galaxies with lower rates of shear, i.e. those with 
rising rotation curves that are closer to solid body ro- 
tation. This model also requires GMCs to be relatively 
long-lived: i.e. several tens of Myr. 

Krumholz & McKee (2005) proposed that the star 
formation rate is regulated by turbulence in a galaxy's 
molecular gas, so that a fixed fraction of the turbulent, 
molecular gas is converted to stars per free fall time of 
the CMC. To recover equation ([2]), one has to assume 
that CMC dynamical timescales are a fixed fraction of 
galactic dynamical timescales and that CMC mass frac- 
tions are constant in different galactic systems (see also 
Wada & Norman 2007). 

1.3. Numerical Studies of GMC Formation and Disk 

SFRs 

The large range of spatial scales involved from GMCs 
to global galactic properties (a range of ~ 10^ from < 10 
parsec-scale clouds to the galactic scale) has made it ex- 
tremely difficult for simulations to encompass GMC for- 
mation on global scales. Work that does perform this is 
either limited to two dimensions (Shctty & Ostriker 2008) 
or has to assume a fixed two-phase medium for the ISM 
(Dobbs 2008). While the clouds do largely reside in the 
plane of the disk (which is in essence a two-dimensional 
system), cloud collisions and feedback can eject gas from 
the surface, a process that controls the pressure of the 
ISM from which the clouds are forming (McKee & Os- 
triker 1977; Cox 2005). Similarly, a fixed phase ISM 
results in clouds being created from gas with a predeter- 
mined structure, from which is it harder to discern the 
main physical effects controlling their formation and evo- 
lution. Three-dimensional models that contained a self- 
consistent multiphase atomic ISM on global scales were 
performed by Tasker & Bryan (2006, 2008) and Wada 
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& Norman (2007), but either did not reach the resolu- 
tion needed to probe below the most massive GMCs (e.g. 
the simulations of Tasker & Bryan (2008) had a limit- 
ing resolution of 25-50 pc), or only considered the inner 
galaxy (e.g. the simulations of Wada & Norman (2007) 
were restricted to r < 2.56 kpc). Robertson & Kravtsov 
(2008) presented a global SPH simulation, with similar 
effective resolution to the simulation of Tasker & Bryan 
(2008) (Kravtsov, private comm.) of star formation in 
disk galaxies including a multiphase ISM, and subgrid 
models for various forms of stellar feedback and molecu- 
lar hydrogen formation and destruction (see also Gnedin, 
Tassis, & Kravtsov 2008). 

Local models can resolve down to smaller scales and, 
if set up in a shearing box, can approximate the effects 
of galactic shear on cloud formation and evolution (e.g. 
Kim & Ostriker 2001, 2006, 2007). Other examples of 
local models have focussed on cloud formation from im- 
posed colliding flows (e.g. Heitsch et al. 2008), including 
the effects of magnetic fields (Heitsch, Stone, & Hart- 
mann 2009). Other groups have studied cloud formation 
in the context of the local interplay between supernova 
feedback and a turbulent, multiphase ISM (e.g. Slyz et 
al. 2005). Simulations including the nonequilibrium for- 
mation and destruction of H2 molecules have been carried 
out by Glover & Mac Low (2007a,b). 

While these local models are very important for study- 
ing many aspects of GMC formation and evolution, they 
cannot explore the global evolution of GMCs as they or- 
bit through the disk or measure how GMC (and star 
formation) properties are related to the more readily- 
observed global galactic properties such as mean gas 
mass surface density. Trends with galactocentric radius 
are also easier to study in global models. 

In this paper we present results from a three- 
dimensional global galaxy (^ 32 kpc box containing 
a gravitationally unstable galactic disk with diameter 
20 kpc) simulation that tracks the formation and evo- 
lution of clouds in a self-consistent multiphase atomic 
ISM. Our main simulation has a limiting resolution of 
7.8 pc. As discussed below, we do not at this stage in- 
clude the physics of star formation and feedback. We 
also do not include magnetic fields, cooling below 300 K, 
or any treatment of the formation and destruction of 
molecules. Our definition of "GMCs", which we also 
refer to interchangeably as "clouds" , is discussed in de- 
tail below, but basically involves a density threshold of 
nn > 100 cm~'^. In real galaxies, the structures that cor- 
respond to this definition will typically be composed of a 
mixture of atomic and molecular gas, but we argue that 
the bulk dynamical properties of the clouds can still be 
reasonably well captured by our simulations. 

2. NUMERICAL TECHNIQUES 

2.1. The Code 

The simulations performed in this paper were run us- 
ing Enzo; a three-dimensional adaptive mesh refinement 
(AMR) hydrodynamics code (Bryan & Norman 1997; 
Bryan 1999; O'Shea et al. 2004). Enzo has been used 
previously in galactic disk simulations (e.g. Tasker & 
Bryan 2006, 2008), where it successfully produced a self- 
consistent multiphase atomic interstellar medium (ISM), 
consisting of a wide range of densities and temperatures. 
Grid codes are particularly adept at modeling multiphase 



gases since the grid cells form natural boundaries allow- 
ing the gas to evolve with a range of temperatures, densi- 
ties and pressures. Particle-based techniques also strug- 
gle to resolve fluid instabilities (Tasker et al. 2008) and 
can suffer from over-mixing problems unless specific al- 
gorithmic steps are taken. 

Enzo evolved the gas using a 3D version of the Zeus 
hydrodynamics algorithm (Stone & Norman 1992). This 
scheme uses an artificial viscosity term to model shocks 
and the variable associated with this, the quadratic arti- 
ficial viscosity, was set to 2.0 (the default) for all simula- 
tions. Radiative gas cooling followed the solar metallicity 
cooling curve of (Sarazin & White 1987) down to temper- 
atures of r = 10"' K and extends down to T = 300 K us- 
ing rates from (Rosen & Bregman 1995). These tempera- 
tures take us to the upper end of the atomic cold neutral 
medium (Wolfire et al. 2003). In this study we do not in- 
clude the formation and destruction of molecules or any 
cooling processes below 300 K. By a combination of dust 
and molecular cooling, gas in real GMCs reaches tem- 
peratures of ~ 10 K, more than an order of magnitude 
below our minimum radiative cooling temperature. How- 
ever, since we are not able to resolve the detailed internal 
structure of clouds and their internal turbulence (a typ- 
ical cloud with a diameter of ~ 100 pc would only have 

13 cells across each linear dimension in our highest res- 
olution simulation) and since we are not including mag- 
netic pressure support, we view our temperature floor of 
300 K as being equivalent to imposing a minimum ID 
signal speed equal to the sound speed Cg = {'jP/pY^^ = 
(7fcr/^)i/2 = 1.80(r/300K)i/2kms-i, where 7 = 5/3, 
/i — 1.273mp (for an assumed tiho — O.lnn)- This sig- 
nal speed is somewhat smaller than observed internal 
velocity dispersions of GMCs. In fact we will see that, 
even though our simulation does not include feedback, 
our simulated GMCs typically attain internal nonther- 
mal (i.e. turbulent) velocity dispersions much greater 
than this minimum. 

The galaxy was modeled in a 3D simulation box of 
side 32 kpc with isolated gravitational boundary condi- 
tions and outflow fluid boundaries. For our main, high- 
resolution simulation, the root grid was 256"^ with an ad- 
ditional four levels of refinement, producing a minimum 
cell size of 7.8 pc. To examine how simulation results 
depend on resolution, we also carried out medium and 
low-resolution runs with a minimum cell sizes of 15.6 pc 
and 31.2 pc, respectively. Refinement of a cell occurred 
when the Jeans length drops below four cell widths, in 
accordance with the criteria suggested by (Truelove et al. 
1997) for resolving gravitational instabilities. Resolution 
of the fragmentation is discussed further in Section 12.41 
The disk was allowed to evolve for 324 Myr, about 1.3 
orbital periods at r = 8.0 kpc (see Section [2?2|) and sub- 
stantially longer than the local free fall time of the initial 
mid-plane gas at this location, which was about 60 Myr. 

In this paper, the first of a series, we examine the for- 
mation and evolution of GMCs in a flat rotation curve 
galactic disk without the presence of star formation and 
feedback mechanisms, which we defer to future papers. 
Since the amount of mass removed by star formation 
is expected to be only a few percent per local free fall 
time (Zuckerman & Evans 1974; Krumholz & Tan 2007), 
the gas depletion over the course of the simulation time 
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would be relatively modest. The structure of the ISM 
and the GMC population we derive can be regarded as 
that resulting from the limiting case when there is very 
little coupling of stellar feedback (including FUV heat- 
ing) to the ISM. Gravitational scattering of bound clouds 
in a shearing disk and dissipative cloud collisions are the 
processes that will act in our simulations to extract or- 
bital energy and regulate the ISM structure. 

2.2. The Initial Structure of the Galactic Disk 

To mimic the present-day Milky Way, the simulated 
galaxy was set up as an isolated disk of gas orbiting in a 
static background potential which represented both dark 
matter and a stellar disk component. This minimized ef- 
fects from the evolution of the galaxy on the interstellar 
structure, allowing us to investigate the formation and 
evolution of GMCs in a steady-state environment. The 
background potential was chosen to produce a flat rota- 
tion curve at r ^ Tc, where the core radius was set to 
be Tc = 0.5 kpc. The form of the potential is (Binney & 
Tremaine 1987): 



^ = l^lo In 



1 



(9) 



where Vcfi is the constant circular velocity in the limit 
of large radii, here set equal to 200 km s^^, r and z are 
the radial and vertical coordinates respectively, and the 
axial ratio of the potential field is q,j, = 0.7. The form of 
the circular velocity of the disk, Vc, is then given by: 



(10) 



The initial radial profile of the gas mass surface den- 
sity in the disk was chosen to give a constant value 
of the Toomre Q parameter for gravitational instability 
(Toomre 1964): 

■^'^ (11) 



where ag, the ID velocity dispersion of the gas, is equiv- 
alent to the sound speed Cs for the case of an razor thin 
disk with only thermal pressure. The particular choices 
of Q at different radii are discussed below. For our sim- 
ulated gas disk we define ag = (a^^ + c^)^/^, averaged 
over the mass in particular regions of the disk, where Cnt 
is the ID velocity dispersion of the gas motions in the 
plane of the disk after subtraction of the circular veloc- 
ity, i.e. representing nonthermal motions. Our disks are 
initialized with ant — 0. 

The initial vertical profile of the gas was set pro- 
portional to sech^{z/zh), where Zh is the vertical scale 
height, which was assumed to vary, i.e. increase, with 
galactocentric radius based on observations of the HI in 
the Milky Way presented in (Binney & Merrifield 1998). 
At the solar radius of 8 kpc, Zh = 290 pc. Then we have 
Sff = J^^Posech^{z/zh)dz = 2poZh, where po is the 
midplane (z = 0) density, so that the gas distribution 
becomes: 



(12) 




tg actual cloud position 

\.J ^1 predicted cloud position 
t.^ actual cloud position 



The initial disk profile is divided radially into three 
sections. In our main region of interest, between radii 



Fig. 1. — Diagram to illustrate how clouds are tracked from 
one output timestep to the next, typically 1 Myr later. Panel (A) 
shows a simple example: the predicted position (center of mass) of 
cloud 1 at time ti is shown by the black dotted circle. A volume 
of radius 50 pc is searched for clouds present at ti and the nearest 
cloud is tagged as being cloud 1 at time ti. Panel (B) shows a 
more complex scenario where two clouds present at time to merge 
to form a single cloud at time fi. (diagram is illustrative and not 
to scale). 



of r = 2 — 10 kpc (that is, encompassing the part of the 
Galaxy inside the solar circle with radius of 8 kpc) Sg is 
set so that Q = 1 if ct^ = were equal to 6 km s~^, simi- 
lar to the observed velocity dispersion of the ISM (Kenni- 
cutt 1998; Stark & Brand 1989). The actual initial veloc- 
ity dispersion (i.e. sound speed) is set to Cg = 9.0km s~^, 
so that Tinit ~ 7450 K and Q = 1.5. The threshold of 
gravitational instability is crossed by allowing the gas to 
cool. The other regions of the galaxy, from to 2 kpc 
and from 10 to 12 kpc, are initialized in a similar way, 
but are designed to be gravitationally stable with Q — 20 
for a flat rotation curve if Cs were to equal 6 km s~^ (al- 
though note the rotation curve is not flat in the center). 
Their initial temperatures are set to the same value as 
in the main disk region. Even after cooling to the tem- 
perature floor of 300 K, for which Cs = 1.8km s~^, these 
regions remain gravitationally stable, except for a small 
amount of gas that collects at the galaxy center. Be- 
yond 12 kpc, the disk is surrounded by a static, very low 
density medium, that has negligible influence. Only the 
main disk region from 2 to 10 kpc is analyzed in this 
paper, and in fact, due to boundary effects, we typically 
restrict analysis to radii from 2.5 to 8.5 kpc. 

In total, the gas mass was 6 x 10® Mq. Note, this is a 
factor of ten smaller than the disks presented in Tasker 
& Bryan (2008). 

2.3. Defining and Analyzing Giant Molecular Clouds 

The clouds in the galactic disk were located using a 
number density of H nuclei threshold of nn.c = 100 cm~^, 
similar to the mean (volume-averaged) densities of typ- 
ical Galactic GMCs. Also recall that the formation of 
molecules is not being followed in this simulation. When 
we refer to 'clouds' or 'GMCs' wc are describing the gas 
that has achieved densities of njj > ^^H.c in the atomic 
phase. We expect that most of the gas above this den- 
sity would form molecules via surface reactions on dust 
grains, although there could still be substantial ('--' equal 
mass) atomic components present, as are observed in and 
around Galactic GMCs (Wannier et al. 1991; Bhtz 1990). 
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Fig. 2. — Example of cloud tracking in the medium resolution (15.6 pc) simulation. The left panel shows a 2 kpc square patch of the 
galactic plane, 5 kpc from the galactic center, at t = 160 Myr, 15.6 pc thick in the z-direction, with the grey scale showing column density 
through this slice, i.e. proportional to njj in the cells. Contours show the cloud definition threshold density of riH.c = 100 cm~^, while the 
numbered squares and triangles show the center-of-mass position of the clouds located inside the volume of the slice with prograde and 
retrograde rotation, respectively. The right panel shows the same patch at t = 161 Myr. Galactic orbital motion has carried clouds and 
the ISM structure upwards and to the left. The cloud tracking algorithm successfully tracks clouds as they move, and records a merger of 
clouds 9282 and 13550 on the right side of the images. Note the center of mass of cloud 13155 enters the slice from above during this time 
interval, just lower-left of center of the image. 



The process of locating and tracking the clouds in the 
disk at a given time in the simulation is outlined below: 

1. Clouds were identified as peaks in the baryon den- 
sity field with riH > t^-h.c — 100 cm~^. If two peaks 
were < 4 minimum cell widths apart, i.e. about 
30 pc in our high resolution run, then only the 
higher density peak was retained for cloud defini- 
tion. 

2. Non-peak cells with riu > ?ih.c were assigned to 
the nearest peak that was connected to them by 
cells with riH > "-h.c, i-e. clouds are continuous 
structures. 

Once the cloud had been defined, properties including 
the cloud's center-of-mass, velocity and angular momen- 
tum were calculated. Multiple clouds could exist in the 
same continuous density structure if it contained more 
than one well-separated peak. 

To follow the evolution of the clouds, simulation out- 
puts were analyzed every 1 Myr and the clouds mapped 
between outputs with a tag number assigned to each 
cloud to follow its life through the simulation. To map 
a cloud between times to and ti, the code performed the 
following steps: 

1. Assuming linear motion, a predicted position of 
each cloud's center-of-mass found at time is cal- 
culated for the cloud at ii. 

2. A volume of radius 50 pc, larger than the expected 
deviations from linear motion due to typical accel- 
erations (e.g. over 1 Myr these are about 20 pc 
at r — 2 kpc due accelerations in the galactic po- 
tential) , centered on the predicted position of each 



cloud is searched for clouds present at time ti. If 
multiple clouds are found in this region at ti, the 
nearest one is chosen to be associated with the 
cloud from to . In cases where two or more clouds 
are associated with the same ti cloud, the nearest 
one is matched and the volume around the pre- 
dicted positions of the other to clouds is searched 
for alternative candidates. 

3. If no clouds at ti are associated with a to cloud, 
then a volume with radius equal to 3x the average 
radius of the to cloud is searched. This allows for 
large, extended clouds whose radius may be > 50 pc 
and whose centers may have shifted due to external 
perturbations. 

4. Any to clouds remaining unassigned may have 
merged with neighboring to clouds. A volume of 
radius 2x the average radius of each ti cloud is 
searched for unassigned to clouds. This value was 
chosen to be lower than for the previous step since 
a recently merged cloud is likely to be fairly ex- 
tended. If found, a merger is declared between 
the to cloud previously associated with the ti cloud 
and the unassigned to cloud. The ti cloud inherits 
the tag number from the more massive of the two 
merged clouds. If there were multiple ti clouds 
close to an unassigned to cloud, the closest was 
chosen. Multiple mergers involving more than two 
clouds were possible, though were typically rare. 

5. Any to clouds remaining unassigned after these 
steps are declared to have been destroyed by non- 
merger processes. 

Note that this method assumes that any cloud that 
has been destroyed in close proximity to another cloud 
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has suffered a merger. This assumption works weU in 
the current simulation, which has no stellar feedback, 
but it remains to be tested once feedback processes are 
operating. 

Figure [1] illustrates two possible examples of the cloud 
mapping. The blue circles show clouds at i = to a-nd the 
red circles show clouds located IMyr later at t = ti. In 
panel (A), the predicted position of cloud 1 at t = is 
marked by the black dotted open circle. The surrounding 
region with radius 50 pc is searched for clouds present at 
ti and two are found. The closest of these is identified 
to be cloud 1 at ti whereas the other cloud is either a 
newly formed cloud or a different cloud present at to. 
Panel (B) shows how the algorithm deals with a more 
complex situation involving a merger. Here, a cloud at 
ti is associated with both clouds 2 and 3. The code 
checks there are no unassigned clouds at ti that are a vi- 
able match for cloud 2 or 3 but finding none, it declares 
this a merger event. Note, that the cloud detection algo- 
rithm does not allow distinct clouds to be found within 
four minimum cell widths of each other. 

Figure [2] shows how the cloud identification and track- 
ing work in the simulation. The two images show the gas 
density in a one-cell thick slice of a square patch of the 
galaxy disk midplane 2 kpc across and 5 kpc from the 
galactic center, taken IMyr apart. Contour lines mark 
the density threshold of nn.c > 100 cm~'^ and the squares 
and triangles (prograde and retrograde rotators, see i j4.2p 
show the center-of-mass of each of the clouds identified, 
with their tag number written below. Areas where con- 
tours or a peak appear to be visible without a cloud are 
due to the cloud's center-of-mass being above or below 
the slice shown. The cloud tracking algorithm correctly 
maps clouds between the two outputs and shows a merger 
of clouds 9282 and 13580 on the right-hand side. 
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Fig. 3. — Mean density of cells as a function of their Jeans 
number, Nj = Ax/Xj for low, medium and high resolution runs 
at t = 250 Myr. The results at earlier times are similar. Note most 
of the gas with nji > lOcm"'^ has cooled to the temperature floor 

1/2 

of 300 K, so for fixed cell size Ax, we have Nj oc . This graph 
shows that in the high resolution run fragmentation of structures 
forming out of material with densities up to riu — TOcm""^ is well- 
resolved, and this applies to most clouds forming from ambient 
galactic gas. The fragmentation of the clouds (which have typical 
volume averaged densities of a few hundred) themselves is less well- 
resolved. The effect of resolutio n on the cloud mass function and 
other properties is examined in i|4.2l 



2.4. Resolving Gravitational Collapse 

Finite computational resources impose a limit on the 
number of levels of refinement we are able to employ in 
the simulation. Cells at the highest level of refinement 
are themselves not refined further, which means they can 
grow in mass and density to the point that the Truelove 
et al. (1997) condition to avoid artificial fragmentation 
is no longer satisfied, i.e. that the Jeans length, Aj, as- 
sociated with gas in each cell should be at least 4 times 
larger than the size of the cell. Ax. In Figure [3] we show 
the mean density of cells as a function of their Jeans 
number, Nj — Ax/Xj. This figure indicates that for the 
high resolution simulation the Truelove et al. criterion is 
satisfied for densities up to nn — 70cm~^, which is close 
to the threshold density used to define our GMCs. Thus 
we expect that the gravitational fragmentation process 
that leads to the formation of GMCs is reasonably well 
resolved in our simulation. It is not the goal of this pa- 
per to follow the fragmentation inside GMCs, and for 
this reason a minimum effective sound speed of about 
1.8 km s~^ is imposed (see above). It is possible that 
large GMCs may suffer spurious fragmentation in the 
simulation, however it should be noted that these GMCs 
typically have internal velocity dispersions that are much 
larger than the sound speed (see We examine the 

effect of simulation resolution on the CMC mass function 
and other properties in m.2[ 



3. GLOBAL ISM STRUCTURE AND EVOLUTION 

The evolution of the disk is shown in Figures [J] & [51 
Starting with a velocity dispersion of ag — Cs — 9 km s~^ , 
(i.e. at a temperature of 7450 K) the disk is gravitation- 
ally stable. The gas cools relatively quickly, lowering the 
sound speed to below 6 km s~^, at which point Q < 1 in 
the main disk region. This cooling occurs most rapidly 
in the denser inner regions. 

The first structure to form is an overdense ring near 
the inner boundary of the main disk at r = 2 kpc. The 
formation of this structure is infiuenced by the boundary 
condition of the inner edge of the main disk: i.e. the 
relative lack of gas inside 2 kpc. The Toomre ring insta- 
bility gathers material radially from a scale about equal 
to the Toomre length. At = iTi^GYig/ (the numerical 
coefficient applies for infinitely thin gas disks), i.e. the 
most unstable scale. Similar ring structures form slightly 
later just outside the inner ring and also at the outer 
edge of the disk. These rings then fragment azimuthally 
into clouds. To avoid structures being influenced by the 
boundaries of the main disk, we restrict our cloud anal- 
ysis to clouds that form between 2.5-8.5 kpc. 

The fragmentation of the main disk begins at the inner 
part of this region, not only because the cooling occurs 
faster in this denser region, but also because once gravi- 
tationally unstable structures have formed, being denser 
they collapse with shorter free-fall times. By t ~ 140 Myr 




Fig. 4. — Evolution of the galactic disk. Images are 20 kpc across and show the disk gas mass surface density, Eg (integrated vertically 
over \z\ < 1 kpc) at f = 50, 100, 200 and 300 Myr. The formation of rings via the Toomre instability is evident at earlier times. These rings 
fragment into individual clouds, which then suffer interactions via galactic differential rotation. The properties of the clouds in this fully 
fragmented stage (t > 140 Myr) are the focus of this paper. 
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Fig. 5. — Galactic disk azimuthally-averaged (60 pc wide annuli) radial profiles and their evolution; from left to right: (a) gas mass 
surface density, Eg = Jl^^^^^ p{z)dz, (b) ID gas velocity dispersion, ag, (mass-weighted average over —1 kpc < z < 1 kpc utilizing only 
disk plane velocity components), (c) gas temperature, T, (mass-weighted average over —1 kpc < z < I kpc), (d) Toomre Q parameter, 
evaluated using Eg and erg. 



the main region of the disk, i.e. out to ~ 8.5 kpc, has 
fully fragmented. 

Once formed, clouds start to interact as differential 
rotation in the disk brings them into contact with one 
another, and again this evolution occurs more rapidly in 
the inner regions of the disk. Gammie, Ostriker, & Jog 
(1991) investigated this process analytically and with nu- 
merical integrations of binary collisions and cloud inter- 
actions, finding interaction times that are a fraction of 
the orbital time. We investigate the cloud collision time 
in more detail in 

Outside the main disk region, the low density gas re- 
mains stable. The very center of the disk, r < 2 kpc, does 
show some development of overdense structures. This is 
due to the infall of gas to the galactic center, most likely 
due to numerical viscosity resulting from the use of a 
Cartesian grid to simulate small-scale circular gas mo- 
tions. These inner regions are excluded from our analy- 
sis. 

Azimuthally-averaged radial profiles of Sg, Ug, T, and 
Q are shown in Figure [51 Note that here we show 

Eg — J^^^^^ p{z)dz, which is effectively equal to the full 
vertical mass surface density through the simulation box, 
CTg is a mass- weighted average over —1 kpc < ^ < 1 kpc 
utilizing only disk plane velocity components, T is a 
mass- weighted average over —1 kpc < z < 1 kpc and 



Q makes use of Eg and Ug via equation ([3]) . 

In the main disk, away from the boundaries, the mean 
value of Eg at a particular radius does not change very 
much during the course of the simulation. Small scale 
fluctuations due to ring instabilities are prominent at 
early times, but become smoothed out as the rings frag- 
ment into clouds, which then gravitationally scatter off 
of each other. 

During the early stage of the simulation, the gas cools 
rapidly, which causes Q to drop below the threshold for 
instability. The gas is also contracting vertically towards 
the disk midplane. After about 50 Myr, the bulk prop- 
erties of the main disk remain relatively constant for 
some time, since the temperature floor of 300 K has been 
reached for much of the disk material. At later times, af- 
ter the disk has fragmented into clouds and they start 
interacting gravitationally, there is a signiflcant increase 
in the velocity dispersion of the gas and some heating. 
At late times the Toomre stability parameter rises to 
values above unity, because of the strong gravitational 
scattering of clouds. 

Figure [6] shows the probability distribution function 
(PDF) for gas density. The top panel shows the volume- 
weighted PDF, evaluated over a volume extending radi- 
ally from 2.5 to 8.5 kpc and ±1 kpc above and below 
the disk midplane. The mass-weighted PDF is shown in 
the bottom panel. These flgures show the relatively fast 
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Fig. 6. — Top panel: Evolution of volume-weighted probabil- 
ity distribution function (PDF) for gas density evaluated from 
2.5 kpc < r < 8.5 kpc and —1 kpc < z < 1 kpc. The jagged nature 
of the Myr PDF is due to finite root grid resolution aliasing of 
the initial conditions, and does not influence the later evolution. 
Higher and lower density regions develop due to gravitational frag- 
mentation. Note there is relatively little evolution between 200 
and 300 Myr. The volume fractions of GMCs (i.e. cells with 
riH > 100 cm"3 = 3.46M0 pc-^) at t = 100,200,300 Myr are 
(5.82,4.28,4.12) X lO"**, respectively. The smooth solid line is a 
log-normal fit to the portion of the PDF at cloud densities (see 
text). Bottom panel: Mass- weighted PDF for gas density evalu- 
ated over the same simulation volume. Again, these distributions 
show relatively little evolution between 200 - 300 Myr, and the 
GMC mass fractions at t = 100,200,300 Myr are 0.492, 0.688, 
0.685, respectively. 



evolution from the initial conditions caused by the early 
cooling and fragmentation. Evolution after 100 Myr pro- 
ceeds more slowly: there is very little change from 200 
to 300 Myr. 

Figure [6] also shows a fit of a log-normal distribution, 
p(lna;)dln.T = (27r(TpQp)^^/^exp(— 0.5CTppp[lnx — Inx]^), 
where x = p/p, to the volume- weighted PDF at the den- 
sities relavant to clouds (see also Wada & Norman 2007; 
Tasker & Bryan 2008). Since we are only fitting to a 



portion of the PDF, here we are only interested in the 
width of the distribution, apDF, not the normalization. 
We find apDF = 2.0. Following the empirical relation 
(jppip = ln[l -I- (3Al^/4)] derived from analysis of simula- 
tions of isothermal, non-self-gravitating supersonic tur- 
bulence (Padoan, Nordlund & Jones 1997, Padoan & 
Nordlund 2002, Krumholz & McKee 2005), where M 
is the ID Mach number, we estimate A4 — 8.5. For a 
sound speed of 1.80km s~^, this corresponds to a velocity 
dispersion of 15kms~^, about 50% larger than the typ- 
ical internal velocity dispersions of clouds (! j4.2[) or the 
disk- mass-averaged velocity dispersions (Fig. [5]). This 
moderate discrepancy may be due to self-gravity skewing 
the high-side of the PDF and/or the eff'ects of shearing 
streaming motions in the disk, which are removed from 
the disk-averaged velocity dispersions. 

The density-temperature phase space of the ISM is 
shown in Figure [T) In the top row, the contours are 
related to the volume in the simulation at the given den- 
sities and temperatures. After disk fragmentation, most 
of the volume is at low densities, rifj ~ 10^^ cm^'^, and 
high temperatures, T ~ 10^ K. The temperature fioor 
of the cooling curve at 300 K is evident on the left hand 
side of these diagrams: most of the GMC material is 
at this effective temperature, i.e. has an effective sound 
speed of 1.8 km s~^, and these clouds occupy very little 
volume. Note that cooler temperatures are possible via 
adiabatic cooling. In the bottom row, the contours are 
related to the mass in the simulation at the given den- 
sities and temperatures. Most mass is in high density, 
Tjjj ~ 1 — 1000 cm"'^, structures, including our defined 
"GMCs" with nn > 100 cm'^. 

The typical local Milky Way total diffuse ISM pressure 
is about 2.8 x lO** K cm""^ and its thermal components 
are about an order of magnitude smaller (Boulares & 
Cox 1990). These pressure are shown by straight lines in 
Figure [71 Our simulated diffuse ISM is at significantly 
lower pressures compared to the observed Milky Way 
pressures. This is not surprising since this simulation 
does not include feedback from star formation, including 
FUV heating, stellar winds, ionization and supernovae. 
Nevertheless, much of the volume of the simulated ISM is 
in approximate pressure equilibrium. The pressure is set 
by energy input from hot gas produced in shocks result- 
ing from cloud-cloud collisions. GMCs in the simulation 
are at much higher pressures than the diffuse ISM, due 
to their self-gravity (see below). In fact the thermal pres- 
sure of the cloud threshold density at the minimum cool- 
ing temperature is about equal to the mean total pressure 
in the local Milky Way ISM observed by Boulares & Cox 
(1990). 

4. GMC PROPERTIES AND EVOLUTION 

4.1. GMC Formation and Merger Rates 

As the disk becomes gravitationally unstable, local 
patches increase in density, many of which reach our 
threshold value of nH,c = 100 cm~'^, at which point 
we recognize them as a GMC. The formation rate of 
these clouds over the course of the simulation is shown 
in Figure [3 There is a burst of cloud formation asso- 
ciated with the initial fragmentation of the disk, with 
the rate peaking at about 400 Myr~^ at a simulation 
time of i = 100 Myr. By t = 140 Myr, cloud forma- 
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Fig. 7. — Density versus temperature contour plots showing the distribution and evolution (100, 200, 300 Myr from left to right) of gas 
volume (top row) and mass (bottom row) in the galaxy disk for 2.5 kpc < r < 8.5 kpc and —1 kpc < 2 < 1 kpc. The temperature floor 
of 300 K is evident: most "GMC" gas has this effective temperature, corresponding Cg ~ 1.8kms~^. Lower temperatures are possible via 
adiabatic cooling. Most of the volume of the ISM exists at lower densities and higher temperatures, with different phases in approximate 
pressure equilibrium. In each panel, solid lines show estimates of the total pressure in the Milky Way, Ptot/k = 2.8 X 10* K cm""^ (top 
line), the total thermal pressure, Ptl,/fc = 0.36 X 10" K cm-3 (middle line), and the thermal pressure excluding the hot gas component, 
Ah.nohot/fc = 0.14 X lO'* K cm-3 (bottom line) (Boulares & Cox 1990). Since our simulations do not yet include feedback (e.g., FUV 
heating, supernovae) or nonthermal pressure components (e.g., magnetic fields, cosmic rays), it is not surprising that our diffuse ISM 
pressure is 1-2 orders of magnitude smaller than the observed Milky Way value. 



tion settles down to a slower and nearly constant rate of 
~ 200 Myr^^. This is partly due to there being a smaller 
non-cloud diffuse gas reservoir and partly due to heat- 
ing of the disk by gravitational interactions of the cloud 
population. We are most interested in the cloud and ISM 
properties after t = 140 Myr, since before this time their 
properties are influenced by the artificially smooth initial 
conditions and their fragmentation via the Toonire ring 
instability. 

Figure[8]also shows the evolution of the total number of 
clouds and the merger rate of clouds. This information is 
used to calculate the mean merger time of clouds, ^merger, 
relative to their orbital time, as a function of galactocen- 
tric radius (Figure (HI). This calculation assumes mergers 
are all binary mergers (the non-binary merger fraction 
is indeed negligibly small). From Figure [9] we see that 
the average merger time settles to be a small fraction, 
~ 0.2 of the orbital time, with only a modest depen- 
dence on galactocentric radius and simulation time for 
t > 140 Myr. This confirms the results of Tan (2000), 
who estimated tmcrgcr/iorbit — 0.2 based on the simpli- 
fied calculation of cloud orbits and binary interactions 
by Gammie et al. (1991). Note that the orbital time for 



Vc = Vc.Q = 200 km s ^ is 

iorbit = 123(r/4 kpc)(i;c/200 km s'^)-^ Myr, (13) 

so that the typical merger time is only ~ 25 Myr, which 
is relatively short compared to traditional estimates of 
cloud collision time scales of hundreds of Myrs that do 
not allow for the essentially 2D geometry of the GMC 
population, gravitational focusing and that cloud inter- 
action velocities are set by galactic shear and are larger 
(by about a factor of two) than the local cloud ve- 
locity dispersion (see discussion in McKee & Ostriker 
2007). Our derived collision timescales are also shorter 
than many estimates of GMC lifetimes (e.g. McKee & 
Williams 1997; Matzner 2002), which suggests that col- 
lisions are likely to be important even when feedback 
mechanisms are taken into account: i.e. an individual 
GMC is just as likely to have its properties dramatically 
altered by a merger than by a destructive mechanism 
such as supernova or ionization feedback. 

The nearly constant ratio of tmcrgcr/iorbit means that 
the cloud collision rate is tied to the global galactic dy- 
namical timescale and, if cloud collisions are the trigger 
for the majority of star formation, then this provides a 
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Fig. 8.— Evolution of the number of GMCs, TVc, the GMC 
formation rate, dN^ f^-^^/dt, and the GMC merger rate, which is 
essentially equal in magnitude to the rate at which clouds are de- 
stroyed via merging, dNc,incTgc/dt, since virtually all mergers are 
binary mergers. The cloud formation rate shows a burst at about 
100 Myr associated with the initial fragmentation of the disk, which 
then decreases by a factor of about two by 140 Myr, after which it 
remains almost constant for the remainder of the simulation. Note 
that clouds can also be destroyed by disruptive collisions and by 
leaving the (inner) boundary of the analyzed volume. Thus in the 
fully-fragmented phase {t > 140 Myr) the overall cloud destruction 
rate is greater than the formation rate by about 10 clouds per Myr. 
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Fig. 9. — Cloud merger timescales (averaged over 50 Myr inter- 
vals of simulation time) compared to orbital timescale as a func- 
tion of galactocentric radius. The merger rate is driven by differ- 
ential rotation in the disk, and the cross-section grows as clouds 
become more gravitationally bound. After the initial fragmenta- 
tion stage at 100 Myr, the cloud merger timescale settles to small 
values ~ 0.24orbit) with only modest dependence on galactocentric 
radius. The dashed line shows the merger times of clouds with 
Mc > 10^ Mq (i.e. the average time for them to undergo a merger 
with a cloud of any mass, though typically with Mc > 10^ Mq [see 
Fig.lTTkl') evaluated over the interval t = 225 - 275 Myr. 



natural mechanism to explain the global SFR-gas content 
correlations observed by Kennicutt (1998) in which the 
overall star formation efficiency per galactic dynamical 
time is small, while having star formation occur mostly 
on small scales at relatively high efhciency in a highly 
clustered mode (Tan 2000). The relation of star forma- 
tion to cloud collisions in global galaxy simulations will 
be investigated in Paper II. 

The following caveats should be considered in the 
above estimates of the merger timescale. The simulated 
merger rate depends on the mass function and gravi- 
tational boundedness of the GMCs. More massive and 
more gravitationally bound GMCs will experience higher 
merger rates, although for typical cloud mass functions 
most of their collisions will be with lower-mass clouds. 
In Figure [9] we also show the merger timescale for clouds 
with Mc > IO^Mq, which is about a factor of two shorter 
than that of the average cloud. We will see below that 
the simulated GMC population mass function and gravi- 
tational boundedness are similar to observed values. The 
total number of clouds is not fully resolved in these simu- 
lations (the most common cloud mass is severalx 10^ M0 
in the high resolution run at late times). Higher reso- 
lution simulations are likely to resolve larger number of 
clouds, which would increase the merger rate, although 
not necessarily the total gas mass that has been com- 
pressed by coUisional shocks. The inclusion of stellar 
feedback would likely lead to a reduced total mass frac- 



tion in clouds and a reduced gravitational boundedness 
of those clouds that do form, and these effects would 
reduce the merger rate. 

4.2. GMC Properties with Simulation Time 

The cloud mass function is shown in Figure [TUk . Since 
there is no star formation feedback in these models, it 
is difficult to destroy gravitationally bound gas clouds, 
except via their merger into more massive clouds or 
via disruptive collisions. Unbound and weakly bound 
clouds can also be destroyed by shearing tidal forces 
due to grazing and close cloud interactions. As a re- 
sult, we expect the mean cloud mass to grow over 
the course of the simulation, and this behavior is seen 
in Figure [TUk . Note, however, that the peak of the 
cloud mass function is at 6 x 10^ Mq (for uniform 
binning in logMc) and this does not vary significantly 
from t =100 to 300 Myr. Note the initial Toomre 
mass in a Q — 1 disk is Mt = SgA|, — 8.9 x 
106(crg/6kms~i)3(r/4kpc)(uc/200kms-i)^iAfQ, about 
an order of magnitude larger than the peak of the 
simulated GMC mass function. The Jeans mass 
in the disk is Mj = irY^gX^ap (G^^g) 2.2 x 

106(crg/6)kms-i)3(r/4kpc)(z;c/200kms-i)-i(Q/l) Mq, 
where Aj = (Tg/(GSg) is the shortest wavelength permit- 
ting gravitational instability in a thin, nonrotating disk 
(Kim & Ostriker 2001). The gas tends to cool to the 
temperature floor of 300 K, for which <Tg ~ 2 km s~^ and 
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Fig. 10. — Normalized distributions of GMC properties at 100 (dashed lines), 200 (dot-dashed lines) and 300 (dotted lines) Myr after 
the start of the simulation, for which the total number of clouds are 4300, 2770, 1840, respectively. Top left: (a) Cloud mass, Mc- The 
total number of clouds with Mc > 10^ Mq are 625, 856, 694 for t = 100, 200, SOOMyr, respectively. Also shown are example observational 
results for the power law slopes of GMC mass functions in the Milky Way (Williams &: McKee 1997) and M33 (Rosolowsky et al. 2003). 
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indicates «vir = 1- Middle right: (e) Vertical (2) component of specific angular momentum, jz. The vertical line indicates the value of 
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thus Q ~ 1/3, resulting in Mj ~ 3 x 10"* Mq. The ef- 
fect of simulation resolution on the cloud mass function 
at t = 250 Myr is shown in Figure [Tib. The number of 
GMCs appears reasonably well-resolved down to masses 
of about 2 X 10^ Mq. 

The simulated cloud mass function can be compared 
to the observed mass spectrum of GMCs in the Milky 
Way and M33, which can be fitted by a power law 
dNcl(d\a.Mc) oc M^"'. In the Milky Way, ac is ob- 
served to be ~ 0.6 - 0.8 (Wilhams & McKee 1997), 
whereas a steeper index of ~ 1.6 appears to hold in M33 
(Rosolowsky et al. 2003). After collision and agglom- 
eration processes have had enough time to operate, i.e. 
by i = 200 Myr, the shape of the simulated cloud mass 
function (in the mass range ~ (0.5 — 10) x 10^ Mq rele- 
vant to observations) approaches that seen in the Milky 
Way CMC population, indicating that most of the total 
cloud mass is contained in the more massive clouds. It 
remains to be determined whether these processes would 
still be as important for shaping the cloud mass function 
if stellar feedback processes were also operating, which 
would tend to reduce cloud lifetimes. 

Williams & McKee (1997) also derived the normal- 
ization of the Galactic GMC mass function and noted 
there appears to be a truncation above a molecular mass 
of M„ = 6 X 10^ Mq (see eq. [6|). Allowing for an 
equal amount of atomic gas associated with GMCs (Blitz 
1990), this truncation would occur at 1.2 x lO'^ Mq. We 
see that at late simulation times, the cloud mass function 
has grown beyond this value, though there are very few 
clouds with Mc > l(f Mq (see also Figure EJi) . The 
presence of these very massive clouds is likely due to the 
lack of stellar feedback cloud destruction mechanisms in 
these current simulations, especially the ability of stars to 
destroy their natal clouds. Williams & McKee (1997) es- 
timate there are ~ 100-200 inner Milky Way GMCs with 
> 10*5 Mq, and about 1000 with > lO^M©. If one 
allows for an equal mass of atomic gas associated with 
the observed molecular gas of these GMCs, then in com- 
parison we find about 400 clouds with > 2 x 10^ Mq, 
a factor of 2 to 4 higher than the Williams & McKee 
estimate. We do not adequately resolve the number of 
2 X 10^ Mq clouds (Figure fTTk ) to make a useful com- 
parison at that mass scale. These results may indicate 
that the number of massive GMCs in our simulation is 
a factor of a few larger than in the Milky Way. Indeed 
the mass fraction of gas in GMCs in the simulation is 
about 0.69 (Figure 0, somewhat higher than the mass 
fraction implied by the analysis of Wolfire et al. (2003). 
Lack of feedback processes in the simulation is an obvi- 
ous potential cause of this discrepancy. Another poten- 
tial contributing factor is our somewhat arbitrary choice 
of nH,c = 100 cm^^ for the cloud threshold density. Note 
that M33, being a much smaller galaxy than the Milky 
Way, has fewer massive GMCs. 

Figure [TOb shows the average radius of the clouds de- 
fined as Rc,A = (^c/tt)^/^, where Ac is the projected area 
of a cloud in the y-z plane (i.e. as it would be observed by 
an observer embedded in the plane of the galaxy). The 
most common radius for a cloud is around 20-30 pc (for 
uniform linear binning in Rc,a)- By 300 Myr, agglomera- 
tion processes have created a population of larger clouds. 
It is clear from Figure [TUb & b that the numbers of these 



larger, more massive clouds are steadily increasing in 
time, and in this respect a steady state has not been 
reached. Such a steady state likely requires feedback 
processes. Nevertheless, the typical sizes of the cloud 
population are similar to observed sizes of GMCs in the 
Milky Way and other galaxies, such as M33 (Rosolowsky 
et al. 2003). The effect of simulation resolution on the 
distribution of cloud sizes is shown in Figure [Tib . 

Figure [TUb shows the distribution of mass surface den- 
sity of the clouds, which is defined as Ec = Mc/Ac. The 
peak of the distribution (for uniform binning in logMc) is 
at about 300Mqpc~^, and does not change significantly 
during the simulation. Although the number of more 
massive clouds is increasing, these are also larger, and so 
have similar values of Ec. Most clouds are within a factor 
of 3 of this peak value. This surface density is similar to 
the mean value of ~ 200 Mq pc~^ derived by Solomon et 
al. (1987) from a ^^CO survey of the Galaxy. Heyer et al. 
(2008) have derived smaller values ~ 100 Mq pc~^ based 
on better sampled ^^CO surveys. To compare these ob- 
served values, which are based solely on molecular line 
(CO) emission, with our simulated clouds, one should 
also include their associated atomic gas, which may con- 
tribute at about the factor of two level. We conclude 
that, even without disruptive effects of star formation 
feedback, our simulated GMCs attain values of Sc that 
are similar to observed values. This may indicate that the 
dominate source of turbulent pressure support in GMCs 
is injected via cloud collisions and interactions (i.e. tur- 
bulent converging flows) rather than via star formation 
feedback. Note that Joung & Mac Low (2006) found 
that supernova driven turbulence was insufficient to ex- 
plain the observed ISM velocity dispersions. The effect of 
simulation resolution on the distribution of Ec is shown 
in Figure lllb : the peak of the distribution is relatively 
insensitive to resolution. 

Figure [TOH shows the distribution of GMC virial pa- 
rameters (see eq. [7]), avir = ^ctIRc.a/ [GMc), where CTc 
is the mass-averaged ID velocity dispersion of the cloud, 
i.e. (Tc = {cl + crnt,c)"^^^' where Cnt.c is the ID RMS 
velocity dispersion about the cloud's center of mass ve- 
locity (Bertoldi & McKee 1992). A virial parameter 
of unity implies a spherical, uniform cloud with negli- 
gible surface pressure and magnetic fields is virialized, 
so that its total kinetic energy is half the magnitude of 
the gravitational energy. Surface pressure terms and a 
centrally concentrated density distribution tend to raise 
the value of avir corresponding to virial equilibrium (see 
eq. [7] and the discussion above). We see that the dis- 
tribution of cloud virial parameters peaks at ctvir < 1, 
and decreases slightly as the simulation progresses: the 
peak of the virial parameter distribution at 300 Myr is 
about 0.6. This indicates that the clouds in the simula- 
tion are typically gravitationally bound, perhaps some- 
what more strongly bound than observed GMCs, though 
there is still a significant population of unbound clouds. 
The effect of increasing simulation resolution on the dis- 
tribution of Q!vir is shown in Figure [TTH : the peak of the 
distribution moves to larger values and the fraction of 
unbound clouds grows. 

Figure [TUb shows the distribution of the vertical (z) 
component of specific angular momentum, j^, of the 
clouds. The fraction of retrograde (i.e. jz < 0) clouds 
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grows at late times as clouds have had more opportu- 
nities to suffer mergers and close interactions (see also 
Dobbs 2008). The vertical line indicates the value of 
jz of a spherical (~ 110 pc radius) region of the initial 
conditions at galactocentric radius r = 4 kpc containing 
10^ Mq. (Kim et al. 2003 derive an analytic expression 
for jz oc i^i?cA ^'^^ 2D patches of galactic disks.) We 
see that GMCs typically have much smaller amounts of 
angular momentum compared to similar masses of dif- 
fuse gas that would be spread out over larger scales and 
be more susceptible to galactic differential rotation. The 
effect of simulation resolution on the distribution of jz 
is shown in Figure lllb : at higher resolution it becomes 
more sharply peaked towards small value of jz ■ 

Figure \TUi shows the distribution of cloud center of 
mass vertical positions relative to the disk midplane, i.e. 
where the vertical acceleration due to the galactic po- 
tential is zero. The RMS heights of the cloud popula- 
tion at 100, 200, 300 Myr are 13, 25 and 51 pc, respec- 
tively, growing at late times because of cloud-cloud grav- 
itational scattering. Note Stark & Lee (2005) derived a 
vertical scale height of Milky Way GMCs of < 35 pc. The 
effect of simulation resolution on the vertical distribution 
of clouds is shown in Figure II IF . 

Figure [T^] shows the velocity dispersion versus size re- 
lation for clouds at t = 250 Myr. The effects of using the 
full velocity information to derive a mass-averaged ID 
internal velocity dispersion versus just using radial ve- 
locities at a particular viewing angle are examined. We 
also consider the effect of a mass cut, i.e. restricting the 
analysis only to the most massive clouds. Such effects 
are important in influencing the normalization and slope 
of the line-width size relation. Our simulated cloud pop- 
ulation has a similar, though somewhat higher, velocity 
dispersion versus size relation to observed GMCs, such 
as those observed in M33 by Rosolowsky et al. (2003), in 
the Milky Way by Soloman et al. (1987), or the average 
of a number of extragalactic systems, CTc = o'pc(^/pc)"'' 
with = 0.60 ± 0.10 and dpc = 0.44 ± 0.15 kms^^, 
derived by Bolatto et al. (2008). The larger velocity 
dispersion of the simulated clouds at a given size scale 
may be due to the lack of star formation feedback caus- 
ing them to be more tightly gravitationally bound than 
real clouds. A discrepancy may also result due to differ- 
ences in how the effective radius of clouds is evaluated for 
the simulation clouds and the real clouds. The disper- 
sion of the simulated clouds about the best fit relation is 
about 2.2km s~^ (for the sample in the lower-right panel 
of Figure [T2|) . corresponding to about 40% of Uc. Solo- 
man et al. (1987) found a dispersion of about 30% in 
their sample of Galactic GMCs (see also Heyer & Brunt 
2004). 

Figure [13] shows the distribution of the angle, 9, be- 
tween the cloud angular momentum vector and the galac- 
tic rotation axis at different times in the simulation. At 
early times, in the initial fragmentation phase clouds are 
all born with prograde rotation {0° < 9 < 90°), i.e. in 
the same sense as the galactic rotation. At later times, 
cloud interactions lead to a build up of a retrograde 
(90° <9 < 180°) cloud population, which is about 30% 
of the total. 

To compare to observed systems, one needs to account 
for the viewing angle to the galaxy and the fact that only 



line of sight velocities are used to determine rotation. 
Figure [14] shows these effects by plotting the distribu- 
tion of position angles, 9', of cloud rotation axes derived 
from radial velocities along sight lines at different incli- 
nation angles to the simulated galaxy a,t t — 250 Myr, 
when the retrograde fraction of 6 is 0.28. These results 
indicate that the derived retrograde fraction depends on 
inchnation angle: e.g. for the 0° (edge on) case the retro- 
grade fraction infered from 9' is about twice that as mea- 
sured by the actual angular distribution of 9. The figure 
also shows the observational results of Rosolowsky et al. 
(2003), who find about 40% (18 out of 45) of the GMCs 
in M33 have apparently retrograde rotation. Their re- 
sults are broadly consistent with our simulation results 
(note their slightly larger bin size), but we emphasize 
that we can only make a qualititive comparison at this 
stage, since our model galaxy was not set up to mimic 
the details of M33 (e.g. its rotation curve), nor does it 
yet include the effects of stellar feedback. Nevertheless 
cloud formation in a self-gravitating, shearing disk and 
cloud evolution via cloud-cloud interactions appear to be 
promising mechanisms to help explain the observed an- 
gular distributions of cloud angular momentum vector 
orientations. 

4.3. GMC Properties with Cloud Age 

Our method of tracking clouds from one output time 
to the next allows us to study the evolution of cloud 
properties as a function of cloud age, t' . Figure [151 
analogous to Figure [TO] shows how the distribution of 
cloud masses, radii, mass surface densities, virial pa- 
rameters, vertical component of specific angular mo- 
mentum, and vertical positions changes with cloud age 
(t' = - 1, 9 - 10, 49 - 50, 99 - 100 Myr) for clouds born 
after t = 140 Myr, i.e. in the fully-fragmented phase of 
the simulation. Younger, especially newly-formed, clouds 
tend to have lower masses, slightly smaller radii, lower 
mass surface densities, and larger virial parameters (i.e. 
are less gravitationally bound). However, young clouds 
have very similar distributions of the vertical component 
of specific angular momentum and of vertical position 
above and below the disk midplane. 

Figure [TBI shows the distribution of the angles, 9, be- 
tween cloud angular momentum vectors and the galactic 
rotation axis for different cloud ages for clouds born after 
t = 140 Myr in the fully-fragmented stage. The fraction 
of retrograde clouds is about 25% of the total and nearly 
independent of cloud age, though it has risen to about 
30% for 100 Myr old clouds. At times t > 140 Myr, 
clouds are being born from overdense gas that has al- 
ready been disturbed by cloud interactions. At earlier 
times {t < 140 Myr) we find much smaller retrograde 
fractions, ~ 10%. These results, especially for the later 
simulation times, may depend on the absence of feedback 
in the simulation, which would lead to a higher mass frac- 
tion in lower-density gas. Clouds forming from such gas 
are expected to form with a higher degree of prograde ro- 
tation, since the initial conditions of the clouds are more 
dispersed and thus influenced by galactic shear. The in- 
fluence of feedback on these aspects of cloud formation 
will be investigated in a future paper. 
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Fig. 12. — Velocity dispersion versus size relation for clouds at 4 = 250 Myr. The left panels show (Tc, the mass-averaged ID internal 
velocity dispersion of clouds evaluated by averaging over the full 3D velocities. The right panels show CTc as it would be evaluated for 
clouds observed only with line of sight velocities in a galaxy viewed at 52° inclination angle, similar to our view of M33. The top panels 
include all clouds, while the bottom panels include only clouds with Mc > 10® Mq, i.e. those for which we have a better ability to resolve 
internal dynamics. The solid line shows the power law fits, (Tc = a-pc{Rc a/p^)"" , to the simulated cloud population with power law slopes 
Oct = 0.95, 0.71, 0.51, 0.46 (clockwise from top left) and Cpc = 0.283, 0.378, 0.97, 1.92 km s~^ , respectively. The dashed line shows the result 
of Rosolowsky et al.'s (2003) study of massive (10^ Mq < < 10^ Mq) GMCs in M33, with = 0.45 ± 0.02 and o-pc = 0.72 km s'^. 
The dotted line shows the result of Solomon et al.'s (1987) study of Galactic GMCs, with Oct = 0.5 ± 0.05 and o-pc = 0.72 ± 0.07 km s'^. 



5. CONCLUSIONS 

We have presented simulations of GMC formation and 
evolution in a marginally unstable gas disk of a flat rota- 
tion curve galaxy, capturing scales from ~ 20 kpc down 
to < 10 pc, with a fully multiphase atomic ISM. In this 
initial study, we have focused on the limit of a negligible 
star formation rate from the gas and negligible feedback 
from those stars. We imposed a minimum effective sound 
speed of 1.8 kms~^ to mimic nonthermal forms of pres- 
sure support in the dense gas. We did not explicitly track 
molecule formation on dust grains, but rather defined 
GMCs to be those structures with nn > 100 cm^"^, simi- 
lar to the mean volume-averaged densities of observed 
GMCs. Using adaptive mesh refinement we resolved 
fragmentation of gas up to densities of about this cloud 



threshold density. We developed methods to track clouds 
through the simulation, including mergers. 

In spite of the simplicity of the model, a surprisingly 
large number of observed GMC properties are approxi- 
mately reproduced by the simulated cloud population at 
late simulation times, t > 140 Myr, once the disk is fully 
fragmented. These include the distributions of cloud 
mass, size, mass surface density, virial parameter, angu- 
lar momentum, vertical height in the disk, line-width size 
relation and distribution of angles of angular momentum 
vectors with respect to the galactic rotation axis. Many 
of the ISM and cloud population properties approach a 
quasi steady state at late simulation times, indicating 
an approximate balance between gravitational instabil- 
ity, heating of the cloud velocity dispersion via gravita- 
tional scattering (and perhaps some numerical diffusion) 
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Fig. 13. — Distribution of the angle, 0, between cloud angular momentum vectors and the galactic rotation axis at different times during 
the course of the simulation. The shaded bars indicated retrograde rotation, and this population grows with time as more and more clouds 
experience collisions and close interactions. 
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Fig. 14. — Distribution of the position angle, 8', of cloud angular momentum vectors with respect to the galactic rotation axis as 
determined from radial velocities for t = 250 Myr, for which the retrograde fraction of 9 is 0.28, and for inclination angles of 0° (relevant 
to edge on disks including the Milky Way), 30°, 45°, 52° (relevant to M33), 60° and 90°. Black line shows the observational results from 
(Rosolowsky et al. 2003) for 45 clouds in M33 of which 18 are retrograde (as defined by 6'). 



and dissipation via cloud collisions and mergers. The 
latter processes occur very frequently: a cloud typically 
suffers a merger every 1/5 of an orbital time, i.e. about 
every 25 Myr at r = 4 kpc. GMC mergers and colli- 
sions, which can be viewed as being somewhat equivalent 
to turbulent converging flows in a sclf-gravitating and 
shearing gas disk, are thus efficient at injecting kinetic 



energy into the clouds (extracted from orbital energy) 
and maintaining near virial balance of the clouds. This 
mechanism may be as important as, or even more im- 
portant than, stellar feedback, such as momentum injec- 
tion by supernovae, stellar winds or ionization, at keeping 
GMCs turbulent. A comparison of these processes will 
be studied in a future paper. 
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Fig. 15. — Normalized distributions of GMC properties, as described in Figure [Till but now showing results for different cloud ages: 
0-1 Myr (solid lines), 9-10 Myr (dashed lines), 49-50 Myr (dot-dashed lines), 99-100 Myr (dotted lines), for those clouds born after 
t = 140 Myr, i.e. in the fully-fragmented phase. Top left: (a) Cloud mass, Mc. Middle left: (b) Cloud radius, Rc,a = (Ac/tt)^^^. Bottom 
left: (c) Mass surface density, = Mc/Ac- Top right: (d) Virial parameter, Qvir- Middle right: (e) Vertical (z) component of specific 
angular momentum, jz- Bottom right: (f) Cloud center of mass vertical positions, z. 
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Fig. 16. — Distribution of the angle between cloud angular momentum vectors and the galactic rotation axis for different cloud ages for 
clouds born after t = 140 Myr in the fully-fragmented stage. The shaded bars indicated retrograde rotation. The fraction of retrograde 
clouds is about 25% of the total and nearly independent of cloud age. 
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Cloud collisions may also be efficient triggers for star 
formation, since they cause dense gas that is already 
somewhat prone to star formation to be pressurized by 
the ram pressure of the converging clouds. This could 
push magnetically subcritical clumps inside GMCs to 
become supercritical, thus triggering star cluster forma- 
tion. Our result that the cloud collision time is a small 
and approximately constant fraction of the local orbital 
time, lends support to theory of cloud collision induced 
star formation to explain the global Kennicutt (1998) 
star formation relations, as proposed by Tan (2000). The 
correlation of star formation activity with cloud collisions 
will be investigated in Paper II. 

The complexity of the star formation process means 
that theoretical models and numerical simulations must 
be closely tested against observed systems. An impor- 
tant future goal is to test how GMC properties and star 
formation activity depend on global galactic properties 
such as rotation curve shape and velocity normalization, 
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